function [beta export] = estimate_avgwage(M)

  global omega_coeffs seomega_coeffs;
  omega_coeffs = []; seomega_coeffs = [];

  data.firmid = M(:,1);
  data.year = M(:,2);
  data.lnVA = M(:,3);     % Value added
  data.lnY = M(:,11);     % Output
  data.lnqempl = M(:,4);  % Quality adjusted labor from Mincer
  data.Llnqempl = M(:,12);% lag    
  data.lnL = M(:,10);     % Employment
  data.L = exp(data.lnL);
  data.LlnL = M(:,13);    
  data.lnInp = M(:,5);    % Intermediate useage
  data.LlnInp = M(:,14);  
  data.Inv = M(:,6);      % Investment
  data.lnK = M(:,7);
  data.LlnK = M(:,15);
  data.Ex = M(:,8);       % Export dummy
  data.nace = M(:,9);     % Nace 2 digit code

  firmvars = 16;
  data.lnw = M(:,firmvars+16);     % Average wage
  data.Llnw = M(:,firmvars+17);     % Average wage, t-1
  
  % Which labor variable do we want to use
  data.lnlabor = data.lnL + data.lnw;
  data.Llnlabor = data.LlnL + data.Llnw;
  
  % Which outcome variable to use
  data.y = data.lnVA;
  
  data.Inv(data.Inv<=0) = NaN;

  data.n = length(data.y);
  disp('Number of firm-year obs'); 
  disp(data.n);

  % Simple OLS (used as starting values)

  X = [ones(data.n,1) data.lnlabor data.lnK];                  % Regressors
  b_ols = inv((X'*X))*X'*data.y;                       % Coeff. estimates
  se = sqrt(diag((data.y-X*b_ols)'*(data.y-X*b_ols)/(data.n-size(X,2))*inv(X'*X)));      % Standard error
  clear X;

  % ----------------------------------------------------
  % STAGE I
  % ----------------------------------------------------
  clear i lnK l poly1 poly2 poly3 poly4 X lambda b;
  %i = data.lnInp;          % Ackerberg and LP
  i = log(data.Inv);        % Olley-Pakes
  lnK = data.lnK;
  l = data.lnlabor;
  Dyear = dummyvar(data.year-1995);
  data.Dyear = Dyear(:,2:end);

  % Polynomial in (i,k,m)
  poly1 = [data.Dyear data.Dyear.*repmat(i,1,9) data.Dyear.*repmat(lnK,1,9)];
  poly2 = [i i.^2 i.^3 lnK lnK.^2 lnK.^3 i.*lnK i.^2.*lnK i.*lnK.^2];     
  poly3 = [data.Ex data.Ex.*i data.Ex.*lnK];
  %poly4 = [l l.^2 l.^3 l.*lnK l.^2.*lnK l.*lnK.^2 l.*i l.^2.*i l.*i.^2 data.Ex.*l];   % Ackerberg
  
  X = [ones(data.n,1) data.lnlabor poly1 poly2 poly3]; 

  [b,ci,r,rint,stats]  = regress(data.y,X);
  betaL = b(2);
  
  alpha = .05;
  t = tinv(1-alpha/2,data.n-length(b));
  se_adj1 = (ci(:,2)-ci(:,1)) ./ (2*t); % Standard Error

  lambda = X*b - betaL*data.lnlabor;
  
  % ----------------------------------------------------
  % STAGE II
  % ----------------------------------------------------
 
  gmmf = @(beta) Ackerberg_gmm(beta,lambda,data);


  b0 = b_ols(3);
  options = optimset('Display','iter','TolFun',1e-8,'TolX',1e-8,'MaxFunEvals',5e7,'MaxIter',10000);
  [betaK,fval,exitflag,output,grad,hessian] = fminunc(gmmf,b0,options);       

  beta = [betaL betaK omega_coeffs];
    
  pr = data.y - beta(1)*data.lnlabor - beta(2)*data.lnK;
  export = [data.firmid data.year data.nace pr];


end
